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(N 

Q Abstract 

O The problem of detecting changes in the statistical properties of a stochastic system and time series 

^ arises in various branches of science and engineering. It has a wide spectrum of important apphcations 

f— I ranging from machine monitoring to biomedical signal processing. In all of these apphcations the 

H 

^ observations being monitored undergo a change in distribution in response to a change or anomaly in the 

»^ environment, and the goal is to detect the change as quickly as possibly, subject to false alarm constraints. 

a 

g In this chapter, two formulations of the quickest change detection problem, Bayesian and minimax, are 

introduced, and optimal or asymptotically optimal solutions to these formulations are discussed. Then 
^ some generalizations and extensions of the quickest change detection problem are described. The chapter 

^ is concluded with a discussion of applications and open issues. 

in 
in 

O I. Introduction 

^ The problem of quickest change detection comprises three entities: a stochastic process under obser- 

^ vation, a change point at which the statistical properties of the process undergo a change, and a decision 

• ^ 

^ maker that observes the stochastic process and aims to detect this change in the statistical properties of 

H 

^ the process. A false alarm event happens when the change is declared by the decision maker before the 

change actually occurs. The general objective of the theory of quickest change detection is to design 
algorithms that can be used to detect the change as soon as possible, subject to false alarm constraints. 

The quickest change detection problem has a wide range of important apphcations, including biomedical 
signal and image processing, quality control engineering, financial markets, link failure detection in 
communication networks, intrusion detection in computer networks and security systems, chemical or 
biological warfare agent detection systems (as a protection tool against terrorist attacks), detection of the 
onset of an epidemic, failure detection in manufacturing systems and large machines, target detection 
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f =K (0, 1) f =K (0 .1, 1) 




(a) Stochastic sequence with samples from /o ~ A/'(0, 1) 
before the change (time slot 500), and with samples from 
/i ~ A/'(0.1, 1) after the change. 

Fig. 1: Detecting a change in the r 




(b) Evolution of the classical Shiryaev algorithm when 
applied to the samples given on the left. We see that the 
change is detected approximately at time slot 1000. 

n of a Gaussian random sequence. 



in surveillance systems, econometrics, seismology, navigation, speech segmentation, and the analysis of 



historical texts. See Section VII for a more detailed discussion of the applications and related references. 

To motivate the need for quickest change detection algorithms, in Fig. [la] we plot a sample path of 
a stochastic sequence whose samples are distributed as AA(0, 1) before the change, and distributed as 
AA(0.1, 1) after the change. For illustration, we choose time slot 500 as the change point. As is evident 
from the figure, the change cannot be detected through manual inspection. In Fig. fib) we plot the evolution 



of the Shiryaev statistic (discussed in detail in Section pTi]), computed using the samples of Fig. la As 
seen in Fig. [lb] the value of the Shiryaev statistic stays close to zero before the change point, and grows 
up to one after the change point. The change is detected by using a threshold of 0.8. 



We also see from Fig. lb that it takes around 500 samples to detect the change after it occurs. Can we 
do better than that, at least on an average? Clearly, declaring change before the change point (time slot 
500) will result in zero delay, but it will cause a false alarm. The theory of quickest change detection 
deals with finding algorithms that have provable optimality properties, in the sense of minimizing the 
average detection delay under a false alarm constraint. We will show later that the Shiryaev algorithm, 
employed in Fig. [lb] is optimal for a certain Bayesian model. 

Earliest results on quickest change detection date back to the work of Shewhart |[T|, Q and Page |[3| 
in the context of statistical process/quality control. Here the state of the system is monitored by taking 
a sequence of measurements, and an alarm has to be raised if the measurements indicate a fault in the 
process under observation or if the state is out of control. Shewhart proposed the use of a control chart to 
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detect a change, in which the measurements taken over time are plotted on a chart and an alarm is raised 
the first time the measurements fall outside some pre-specified control limits. In the Shewhart control 
chart procedure, the statistic computed at any given time is a function of only the measurements at that 
time, and not of the measurements taken in the past. This simplifies the algorithm but may result in a 
loss in performance (unacceptable delays when in detecting small changes). In ISj], Page proposed that 
instead of ignoring the past observations, a weighted sum (moving average chart) or a cumulative sum 
(CuSum) of the past statistics (likelihood ratios) can be used in the control chart to detect the change more 
efficiently. It is to be noted that the motivation in the work of Shewhart and Page was to design easily 
implementable schemes with good performance, rather than to design schemes that could be theoretically 
proven to be optimal with respect to a suitably chosen performance criterion. 

Initial theoretical formulations of the quickest change detection problem were for an observation model 
in which, conditioned on the change point, the observations are independent and identically distributed 
(i.i.d.) with some known distribution before the change point, and i.i.d. with some other known distribution 
after the change point. This observation model will be referred to as the i.i.d. case or i.i.d model in this 
article. 

The i.i.d. model was studied by Shiryaev Q, |j5|, under the further assumption that the change point is 
a random variable with a known geometric distribution. Shiryaev obtained an optimal algorithm that 
minimizes the average detection delay over all stopping times that meet a given constraint on the 
probability of false alarm. We refer to Shiryaev's formulation as the Bayesian formulation; details are 



provided in Section III 



When the change point is modeled as non-random but unknown, the probability of false alarm is 
not well defined and therefore false alarms are quantified through the mean time to false alarm when 
the system is operating under the pre-change state, or through its reciprocal, which is called the false 
alarm rate. Furthermore, it is generally not possible to obtain an algorithm that is uniformly efficient 
over all possible values of the change point, and therefore a minimax approach is required. The first 
minimax theory is due to Lorden ||6| in which he proposed a measure of detection delay obtained by 
taking the supremum (over all possible change points) of a worst-case delay over all possible realizations 
of the observations, conditioned on the change point. Lorden showed that the CuSum algorithm of [3] 
is asymptotically optimal according to his minimax criterion for delay, as the mean time to false alarm 
goes to infinity (false alarm rate goes to zero). This result was improved upon by Moustakides f!\ who 
showed that the CuSum algorithm is exactly optimal under Lorden's criterion. An alternative proof of 
the optimality of the CuSum procedure is provided in ||8|. See Section IV for details. 
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PoUak |9] suggested modifying Lorden's minimax criterion by replacing the double maximization of 
Lorden by a single maximization over all possible change points of the detection delay conditioned on the 
change point. He showed that an algorithm called the Shiryaev-Roberts algorithm, one that is obtained 
by taking a limit on Shiryaev's Bayesian solution as the geometric parameter of the change point goes 
to zero, is asymptotically optimal as the false alarm rate goes to zero. It was later shown in |10] that 
even the CuSum algorithm is asymptotically optimum under the PoUak's criterion, as the false alarm rate 
goes to zero. Recently a family of algorithms based on the Shiryaev-Roberts statistic was shown to have 
strong optimality properties as the false alarm rate goes to zero. See |jTT| and Section IV for details. 

For the case where the pre- and post-change observations are not independent conditioned on the 



change point, the quickest change detection problem was studied in the minimax setting by 1 10| and in 
the Bayesian setting by p2|. In both of these works, an asymptotic lower bound on the delay is obtained 



for any stopping rule that meets a given false alarm constraint (on false alarm rate in |10| and on the 
probability of false alarm in [I2j), and an algorithm is proposed that meets the lower bound on the 



detection delay asymptotically. Details are given in Section III-B and Section IV-C 



To summarize, in Sections III and IV we discuss the Bayesian and Minimax versions of the quickest 
change detection problem, where the change has to be detected in a single sequence of random variables. 



and where the pre- and post-change distributions are given. In Section VI we discuss variants and 
generalizations of the classical quickest change detection problem, for which significant progress has been 
made. We consider the cases where the pre- or post-change distributions are not completely specified 
(Section |VI-A[ ), where there is an additional constraint on the cost of observations used in the detection 



process (Section VI-B i, and where the change has to detected using multiple geographically distributed 



sensor nodes (Section VI-C I. In Section VII we provide a brief overview of the applications of quickest 



change detection. We conclude in Section VIII with a discussion of other possible extensions and future 
research directions. 

For a more detailed treatment of some of the topics discussed in this chapter, we refer the reader to 



the books by Poor and Hadjiliadis 1 13| and Chow, Robbins and Siegmund 1 14|, and the upcoming book 



by Tartakovsky, Nikiforov, and Basseville 1 15|. We will restrict our attention in this chapter to detecting 
changes in discrete-time stochastic systems; the continuous time setting is discussed in p3| . 
In Table ^ a glossary of important symbols used in this chapter is provided. 
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TABLE I: Glossary 



Symbol 


Definition/Interpretation 


0(1) 




X ^ 0(1) as c ^ Co, if Ve > 0, 35 > s.t., \x\ < e if \c - cq] < 6 


0(1) 




X = 0(1) as c ^ Co, if 3e > 0, 5 > s.t., \x\ < e if c - co < (5 


g{c) ^ h{c) as c — > 


Co 


limc^co = 1 or g{c) = /i(c)(l + o(l)) as c ^- cq 


{Xk} 




Observation sequence 


Stopping time t on 


{Xk} 


I{T=ri} = or 1 depends only on the values of Xi, . . . , X„ 


Change point r,7 




Time mdex at which distribution of observations changes from Jq to ji 


Pn (E„) 




TX "11'!'* / i1 1 ii' 

Probability measure (expectation) when the change occurs at time n 


Poo (Eoo) 




Probability measure (expectation) when the change does not occur 


cSS sup ^ 




Essential supremum of A, i.e., smallest A such that F(A < il) = 1 


DifiWfo) 




T T^i\'f>ro"pnpf> hptwfpn fi inH HpfinpH IP* ■ I Inn' -^^^^^ \ 

JA. J_/ J^i V Cl^CIH-C UCL W cell ^ ]^ ClllLl J LlCllllCU d j ILL* ^ 1 y ^ ^ J 


(x)+ 




max{x, 0} 


ADD(t) 




ADD(r) = Er=oIP(r = ^) lEn[(r-r)+] 


PFA(r) 




pfa(t) = P(r < r) = Er=o iP(r - < r) 


FAR(t) 




FAR(t) = p \ , 


WADD(r) 




WADD(t) = sup ess sup E„ [(t - n)+\Xi, . . . ,X„_i] 

n>l 


CADD(r) 




CADD(r) = sup E„[r - n|r > n\. 

n>l 



II. Mathematical Preliminaries 

A typical observation process will be denoted by sequence n = 1, 2, . . .}. Before we describe the 
quickest change detection problem, we present some useful definitions and results that summarize the 
required mathematical background. For a detailed treatment of the topics discussed below we recommend 
d, d' 03 and |[T8|. 



A. Martingales 

Definition 1: The random sequence {Xn,n = 1,2,...} is called a martingale if E[X„] is finite for all 

n, and for any ki < k2 < ■ ■ ■ < kn < kn+i, 

'^[Xk„+i\Xki, ■ ■ ■ ,XkJ = (1) 

If the "=" in ([T]| is replaced by "<", then the sequence {Xn} is called a supermartingale, and if the 
"=" is replaced by ">", the sequence is called a submartingale . A martingale is both a supermartingale 
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and a submartingale. 

Some important and useful results regarding martingales are as follows: 

Theorem 2.1 ( [^14j): (Kolmogorov's Inequality) Let n = 1,2,...} be a submartingale. Then 

P ( max Xk > < ^iS^ V 7 > 

\l<k<n J 7 

where = max{0,X„}. 

Kolmogorov's inequality can be considered to be a generalization of Markov's inequality, which is given 
by 

E[X+1 

P > 7) < — V 7 > (2) 

7 

As we will see in the following sections, quickest change detection procedures often involve comparing 
a stochastic sequence to a threshold to make decisions. Martingale inequalities often play a crucial role 
in the design of the threshold so that the procedure meets a false alarm constraint. We now state one of 
the most useful results regarding martingales. 

Theorem 2.2 ([ [16]): (Martingale Convergence Theorem) Let {Xn,n, = 1, 2, . . .} be a martingale (or 
submartingale or supermartingale), such that sup„E[|X„|] < oo. Then, with probability one, the limit 
Xoo = limfc_i.oo Xn exists and is finite. 

B. Stopping Times 

Definition 2: A stopping time with respect to the random sequence {Xn,n = 1,2, . . .} is a random 
variable r with the property that for each n, the event {t = n} £ o'{Xi, . . . , Xn), where cr{Xi, . . . , Xn) 
denotes the sigma-algebra generated by {Xi, . . . ,Xn). Equivalently, the random variable which 
is the indicator of the event {r = n}, is a function of only Xi, . . . , Xn- 

Sometimes the definition of a stopping time r also requires that r be finite almost surely, i.e., that 
P(r < oo) = 1. 

Stopping times are essential to sequential decision making procedures such as quickest change detection 
procedures, since the times at which decisions are made are stopping times with respect to the observation 
sequence. There are two main results concerning stopping times that are of interest. 

Theorem 2.3 ( p4^): (Doob's Optional Stopping Theorem) Let = 1,2, . . .} be a martingale, 

and let r be a stopping time with respect to {X^ n = 1, 2, . . .}. If the following conditions hold: 

1) P(r < oo) = 1. 

2) E[\Xr\] < oo. 

3) E[X„I|,->„}] — ^ as n — ^ oo. 
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then 

E[Xr]=E[Xi]. 

Similarly, if the above conditions hold, and if n = 1, 2, . . .} is a submartingale, then 

and if n = 1,2,...} is a supermartingale, then 

E[X^] < E[Xi]. 

Theorem 2.4 ( [17]): (Wald's Identity) Let {Xmn = 1,2,...} be a sequence of independent and 
identically distributed (i.i.d.) random variables, and let r be a stopping time with respect to = 
1,2,...}. Furthermore, define the sum at time n as 

n 

Sn = 

k=l 

Then, if E[|Xi|] < oo and E[t] < oo, 

E[Sr] =E[Xi]E[r] 

Like martingale inequalities, the optional stopping theorem is useful in the false alarm analysis of quickest 
change detection procedures. Both the optional stopping theorem and Wald's identity also play a key role 
in the delay analysis of quickest change detection procedures. 

C. Renewal and Nonlinear Renewal Theory 

As we will see in subsequent sections, quickest change detection procedures often involve comparing 
a stochastic sequence to a threshold to make decisions. Often the stochastic sequence used in decision- 
making can be expressed as a sum of a random walk and possibly a slowly changing perturbation. To 
obtain accurate estimates of the performance of the detection procedure, one needs to obtain an accurate 
estimate of the distribution of the overshoot of the stochastic sequence when it crosses the decision 
threshold. Under suitable assumptions, and when the decision threshold is large enough, the overshoot 
distribution of the stochastic sequence can be approximated by the overshoot distribution of the random 
walk. It is then of interest to have asymptotic estimates of the overshoot distribution, when a random 
walk crosses a large boundary. 

Consider a sequence of i.i.d. random variables {Yn} (with Y denoting a generic random variable in 
the sequence) and let 

n 

Sn — ^ ^ Yk , 

k=l 
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and 

r = inf{n > 1 : 5'„, > 5}. 

The quantity of interest is the distribution of the overshoot St — b. If {Yn} are i.i.d. positive random 
variables with cumulative distribution function (c.d.f.) F{y), then {Yn} can be viewed as inter-arrival 
times of buses at a stop. The overshoot is then the time to next bus when an observer is waiting for a 
bus at time b. The distribution of the overshoot, and hence also of the time to next bus, as 6 — )• oo is a 
well known result in renewal theory. 

Theorem 2.5 ( /|i 7/ ): If {Y^} are nonarithmetic[^ random variables, and F(Y > 0) = 1, then 

POO 

lim F{Sr -b>y) = (E[y])"i / F{Y > x}dx. 

b~^oo Jy 

Further, if E[y2] < oo, then 

lim E{Sr-b) ' ' 



b^oo ' ' 2E[Y] 

When the {Yn} are i.i.d. but not necessarily non-negative, and ¥,[Y] > 0, then the following concept 
of ladder variables can be used. Let 

r+ = inf{n > 1 : 5„ > 0}. 

Note that if r+ < oo, then Sr^ is a positive random variable. Also, if 

r = inf{n > 1 : 5"^ > 6} < oo 

then the distribution of St — b is the same as the overshoot distribution for the sum of a sequence of 
i.i.d. positive random variables (each with distribution equal to that of St^) crossing the boundary b. 
Therefore, by applying Theorem |2.5[ we have the following result. 
Theorem 2.6 ( [^17j): If {Yn} are nonarithmetic, then 

POO 

lim F{St -b>y) = (E[5^J)-i / F{St^ > x) dx. 

b^oo Jy 

Further, if E[Y'^] < oo, then 

E[52 ] 

lim E{St - 6) - + 



b^oo ' ^ ' 2E[StJ' 
Techniques for computing the required quantities involving the distribution of the ladder height 5^^ in 
Theorem |2.6| can be found in (TTJ. 

'a random variable is aritlimetic if all of it probability mass is on a lattice. Otherwise it is said to non-arithmetic. 



October 23, 2012 



DRAFT 



9 



As mentioned earlier, often the stochastic sequence considered in quickest change detection problem 
can be written as a sum of a random walk and a sequence of small perturbations. Let 

n 
k=l 

and 

T = inf{n >l:Zn> b}. 

Then, 



k=l 



Therefore, assuming that E[t] < cxd, Wald's Identity (see Theorem 2.4 1 implies that 

T 



E[Zr] = E 



+ nVr]- (3) 



k 

k=l 

E[T]E[Y]+E[r]r]. (4) 



Thus, 

E[Zr] - E[r]r] 



E\t] 



E[Y] 

b + E[Zr - 6] - E[r]r] 



E[Y] 

If EItJt] and E[Zr — b] are finite then it is easy to see that 

where ~ is as defined in Table H] 

But if we can characterize the overshoot distribution of {Z„} when it crosses a large threshold then 
we can obtain better approximations for E[t]. Nonlinear renewal theory allows us to obtain distribution 
of the overshoot when {r/„} satisfies some properties. 

Definition 3: {ijn} is a slowly changing sequence if 

n'^ max{|r/i|, ■ ■ . , |7?n|} "^°°> 0, (5) 

l.p. 

and for every e > 0, there exists n* and 6 > such that for all n> n* 

P ( max \rjn+k - ??n| > e ) < e- (6) 

\l<k<n5 J 

If indeed {r?n} is a slowly changing sequence, then the distribution of Z^- — b, as b ^ oo, is equal to the 
asymptotic distribution of the overshoot when the random walk Sn = X]fc=i crosses a large positive 
boundary, as stated in the following result. 
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Theorem 2.7 ( [17]): If {1^} are nonarithmetic and {r]n} is a slowly changing sequence then 

lim ¥{Zr -h<x)= lim ¥{Sr -b<x). 

6— >oo fe— >oo 



Further, if Var(y) < cxd and certain additional conditions ((9.22)-(9.27) in |17|) are satisfied, then 



¥\Sl ] 

where C = ^Eis^ ] ' ^ limit of {rjn\ in distribution. 

III. Bayesian quickest change detection 

As mentioned earlier we will primarily focus on the case where the observation process is a 

discrete time stochastic process, with X„ taking real values, whose distribution changes at some unknown 
change point. In the Bayesian setting it is assumed that the change point is a random variable T taking 
values on the non-negative integers, with vr^ = P{r = n}. Let P„ (correspondingly E„) be the probability 
measure (correspondingly expectation) when the change occurs at time r = n. Then, Pqo and Eqo stand 
for the probability measure and expectation when r = oo, i.e., the change does not occur. At each time 
step a decision is made based on all the information available as to whether to stop and declare a change 
or to continue taking observations. Thus the time at which the change is declared is a stopping time r 



on the sequence {Xn} (see Section II-B i. Define the average detection delay (ADD) and the probability 
of false alarm (PFA), as 

oo 

ADD(r) = E[(r-r)+] =5]^„E„ [(r-r)+] (7) 

n=0 

oo 

PFA(r) = P(T<r) = J]7r„P„(T<r) (8) 

n=0 

Then, the Bayesian quickest change detection problem is to minimize ADD subject to a constraint on 
PFA. Define the class of stopping times that satisfy a constraint a on PFA: 

C^ = {t: PFA(t) < a}. (9) 

Then the Bayesian quickest change detection problem as formulated by Shiryaev is as follows. 

Shiryaev's Problem: For a given a, find a stopping time r G Cq. to minimize ADD(t). (10) 



Under an i.i.d. model for the observations, and a geometric model for the change point T, (10 1 can be 



solved exactly by relating it to a stochastic control problem Q, |[5|. We discuss this i.i.d. model in detail 



in Section III-A When the model is not i.i.d., it is difficult to find algorithms that are exactly optimal. 
However, asymptotically optimal solutions, as a — )• 0, are available in a very general non-i.i.d. setting 
|[T2), as discussed in Section [III-B| 
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A. The Bayesian I.I.D. Setting 

Here it is assumed that conditioned on the change point F, the random variables {X„} are i.i.d. with 
probabiUty density function (p.d.f.) /o before the change point, and i.i.d. with p.d.f. fi after the change 
point. The change point F is modeled as geometric with parameter p, i.e., for < p < 1 

7r„ = P{F = n} =p(l I|„>i|, ^0 = (11) 

where I is the indicator function. The goal is to choose a stopping time r on the observation sequence 



{Xn} to solve (10 1. 



A solution to ([TO]) is provided in Theorem 3.1 below. Let X'^ = {Xi, . . . , X„) denote the observations 
up to time n. Also let 

Pn = P(r < n I X^) (12) 

be the a posteriori probability at time n that the change has taken place given the observation up to time 
n. Using B ayes' rule, pn can be shown to satisfy the recursion 

Pn+l = ^{Xn+l,Pn) (13) 

where 

^Xn+l,Pn) = . \ ^ , . (14) 

PnL[Xn+l) + [l-pn) 

Pn=Pn + 0-- Pn)p, L{Xn+i) = /i (X„+i)//o ) is the likelihood ratio, and po = 0. 

Definition 4: (KuUback-Leibler (K-L) Divergence). The K-L divergence between two p.d.f.'s /i and 
/o is defined as 

D{h\\h)= I fiix) log^dx. 

J Jo (a:) 

Note that L'(/i||/o) > with equality iff fi = /o almost surely. We will assume that 

0<I?(/i||/o) <oo. 



Theorem 3.1 ( 11^): The optimal solution to Bayesian optimization problem of ([TO]l is the Shiryaev 
algorithm/test, which is described by the stopping time: 



Ts = inf {n > 1 : > ^a} 



(15) 



if Aa € (0, 1) can be chosen such that 



PFAfr,) = a. 



(16) 
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Proof: Towards solving ( lOi, we consider a Lagrangian relaxation of this problem that can be solved 
using dynamic programming. 



J* = min(^E[(T-r)+] +A/P(T<r)j (17) 
where Aj is the Lagrange multiplier, > 0. It is shown in Q, Q that under the assumption ([16]), there 



exists a A/ such that the solution to ( [T7| ) is also the solution to ( fTO] ). 

Let ©„ denote the state of the system at time n. After the stopping time r it is assumed that the 
system enters a terminal state T and stays there. For n < r, we have B„ = for n < F, and G„ = 1 
otherwise. Then we can write 

v-i 

ADD(r) = IE 5]l{e„=i} 

.n=0 

Furthermore, let Dn denote the stopping decision variable at time n, i.e., Dn = if < r and Dn = 1 



and PFA(t)=E[I|0^=o}]- 



otherwise. Then the optimization problem in pT) can be written as a minimization of an additive cost 
over time: 



J* = minE 



.n=0 



with 



9n{0,d) 



[lt{e=i}I{d=o} + ^/ I{e=o}I{(i=i}] • 



Using standard arguments [19] it can be seen that this optimization problem can be solved using infinite 
horizon dynamic programming with sufficient statistic (belief state) given by: 



1 I 



(F < n I 



which is the a posteriori probability of (12 1. 



The optimal policy for the problem given in ( fTT] ) can be obtained from the solution to the Bellman 
equation: 

J {pn) = min \f (1 - Pn)l{d„=i} + II{d„=0} [Pn + Aj{Pn)] (18) 



where 



Aj{pn) =nJ{HXn+l,Pn))]- 



It can be shown by using an induction argument that both J and Aj are non-negative concave functions 
on the interval [0, 1], and that J(l) = Aj{l) = 0. Then, it is easy to show that the optimal solution for 
the problem in ( [TT] ) has the following structure: 

Ts = inf{A: >l:pn>A]. 
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0.2 ^ 0.4 0.6 0.3 1 50 100 150 

p„ >■ r 

(a) A plot of the cost curves for A/ — 10, p — 0.01, (b) Typical Evolution of the Shiryaev algorithm. Threshold 
fo ~ AA(0, 1), fi ~ AA(0.75, 1) 0.99 and change point T = 100. 

Fig. 2: Cost curves and typical evolution of the Shiryaev algorithm. 



See Fig. 2a for a plot of A/(l — p„) and pn + Aj{pn) as a function of p„. Fig. 2b shows a typical 
evolution of the optimal Shiryaev algorithm. 

We now discuss some alternative descriptions of the Shiryaev algorithm. Let 

Pn 



and 



A. 



Pn 



We note that A„ is the likelihood ratio of the hypotheses 'Hi : F < n' and 'Hq : F > n' averaged over 
the change point: 

Pn 



P(F < n\X'^ 



> n\Xf) 

(i-p)"nr=i/o(^.) 

^ n n 



(19) 



fc=i 



i=k 



where as defined before L{X.i) = Also, i2„ p is a scaled version of A„: 



..71 n 

^ '^^ fc=i i=fc 



(20) 
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Like pn, Rn^p can also be computed using a recursion: 

Rn+l,p = — ; —L{Xn+l), Ron = 0. 

\ — p 

It is easy to see that A„ and p have one-to-one mappings with the Shiryaev statistic p„. 

Algorithm 1 (Shiryaev Algorithm): The following three stopping times are equivalent and define the 
Shiryaev stopping time. 

Ts = inf{n >l:pn>A} (21) 
Ts = inf{n > 1 : A„ > a} (22) 



Ts = inf{n > 1 : Rn,p > -} (23) 

P 



with a = 



We will later see that defining the Shiryaev algorithm using the statistic A„ ([19]) will be useful in 



Section III-B where we discuss the Bayesian quickest change detection problem in a non-i.i.d. setting. 
Also, defining the Shiryaev algorithm using the statistic i?„ p (20 1 will be useful in Section IV where we 



discuss quickest change detection in a minimax setting. 

B. General Asymptotic Bayesian Theory 

As mentioned earlier, when the observations are not i.i.d. conditioned on the change point V, then 



finding an exact solution to the problem (10 1 is difficult. Fortunately, a Bayesian asymptotic theory can 



be developed for quite general pre- and post- change distributions |[T2|. In this section we discuss the 



results from 1 12 1 and provide a glimpse of the proofs. 

We first state the observation model studied in |12]. When the process evolves in the pre-change 
regime, the conditional density of Xn given is /o,n(X„| X""^). After the change happens, the 

conditional density of Xn given X"^^ is given by 

As in the i.i.d. case, we can define the a posteriori probability of change having taken place before 
time n, given the observation up to time n, i.e., 

<n\Xl) (24) 



with the understanding that the recursion ( [T3] ) is no longer valid, except for the i.i.d. model. 

We note that in the non-i.i.d. case also A„ = -^z^ is the likelihood ratio of the hypotheses "Hi : F < n" 



and "i^o : F > n". If 7r„ = P{F = n}, then following (19i, A„ can be written for a general change point 
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distribution {7r„} as 

An 



(1 -Pn) 

P(r < n|Xf ) 



> n|Xf ) 

>^)nr=i/o.(^ii^r') 



' k=l i=k 



where 



fi,i{Xi\xi 

/o,i {Xi I 

If r is geometrically distributed with parameter p, the above expression reduces to 



^ ^' k=l i=k 

In fact, An can even be computed recursively in this case: 

1 + A, 
1-p 



A„+i = ^r^^Yn+i (25) 



with Ao = 0. 

In ijTll, it is shown that if there exists q such that 



^ n+t 

- ^ Yi ^ q a.s. Pn when t — )- oo Vn (26) 

i=n 

(q = D(/i||/o) for the i.i.d. model), and some additional conditions on the rates of convergence are 
satisfied, then the Shiryaev algorithm ( |2T] ) is asymptotically optimal for Bayesian optimization problem 
of ([To]) as a — 0. In fact, Tj minimizes all moments of the detection delay as well as the moments of 
the delay, conditioned on the change point. The asymptotic optimality proof is based on first finding a 
lower bound on the asymptotic moment of the delay of all the detection procedures in the class Ca, as 



a — )■ 0, and then showing that the Shiryaev stopping time ( |2T] ) achieves that lower bound asymptotically. 

To state the theorem, we need the following definitions. Let q be the limit as specified in (|26]l, and let 
< e < 1. Then define 

f ^ n+t 



{, n+t 
t>l■.\-Y,y^ 
i=n 



Thus, Te*^"^ is the last time that the log likelihood sum X]r=^n ^* ^^^^^ outside an interval of length e 
around q. In general, existence of the limit q in ([26]) only guarantees P„(Te*^"^ < oo) = 1, and not the 
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finiteness of the moments of TfS'^\ Such conditions are needed for existence of moments of detection 
delay of t^. In particular, for some r > 1, we need: 

E„[re(")]'^ < oo for all e > and n > 1, (27) 

and 

oo 

^„E„[T, < oo for all e > 0. (28) 

Now, define 



n=l 



d 



lim 

71— >00 



logP(r > n) 



n 



The parameter d captures the tail parameter of the distribution of T. If F is 'heavy tailed' then d = 0, 
and if F has an 'exponential tail' then d > 0. For example, for the geometric prior with parameter p, 

d= |log(l-p)|. 



Theorem 3.2 ( ply): If the likelihood ratios are such that ([26]) is satisfied then 
1) If a = = ^^f^, then Tj as defined in (2]_i belongs to the set Ca- 



2) For all n > 1, the m-th moment of the conditional delay, conditioned on F = n, satisfies: 

(II I \ 
f 1 + 0(1)) as a ^ 0. 
q + d J \ J 



(29) 



3) For all n > 1, if ( p7| ) is satisfied then for all m < r, 

log a 



E„[(rs - n) 



< 



1 + 0(1) I as a ^ 00 



q + d 

I log a 

q + d 

4) The m-th (unconditional) moment of the delay satisfies 

inf E[(r-F)+1'" > 
r6C„ ' ^ ~ \ q + d 



1 + 0(1) ) as a — > if a = Oc 



I — a 



a 



(30) 



1 + 0(1) as a 0. 



(31) 



5) If ( [27| ) and ( [28] ) are satisfied, then for all m < r 

log a 



E[(rs 



< 



q + d 
I log 



1 + 0(1) ) as a — ;> 00 



1 + 0(1) I as a — )• if a = Oq, 



1 — a 



(32) 



q + d J \ 'J a 

Proof: We provide sketches of the proofs for part 1), 2) and 3). The proofs of 4) and 5) follow by 
averaging the results in 2) and 3) over the prior on the change point. 
I) Note that 
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Thus, Aa = 1 — a would ensure PFA(rs) < a. Since, Cq, = jz^^ we have the result. 
2) Let La be a positive number. By Chebyshev inequality. 



Kn[(r-n)- 



P„((r-nr>L^)< 
This gives a lower bound on the detection delay 

E„[(r-n)+]- > L-P„((r-nr>L-) 
= L^FniT-n> La). 

Minimizing over the family Ca, we get 



inf En[{T-n)+r>L 

tGCo 



inf FJt -n> Lc 

reCo, 



Thus, if 



inf P„(t - n > L^) — )• 1 as a — )• 



(33) 



then LJJ* is a lower found for the detection delay of the family Ca- It is shown in |12| that if (26 1 



is satisfied then ( [33] ) is true for Lq = (1 — e) ^ '°f ?^ for all e > 0. 



3) We only summarize the technique used to obtain the upper bound. Let {Sn} be any stochastic 
process such that 

Sn, 

g as n — 7- oo. 



Sn 

n 



Let 



and for e > 



r] = inf{n > I Sn > b} 



Te = sup{n > 1 : 



b 



n 



>6}. 



First note that Srj^i < b < Srj. Also, on the set {k > T^}, \Sn/n — q\ < e for all n > k. The event 
{\Sn/n — q\ < e} implies n < Using these observations we have 



rj — 1 



= iV- l)I{r,-l>TJ + (f? - l)I{r,-l<TJ 
< iV-l)hr,~l>T.}+T, 

b 



< 



If E[Te] < oo, and because e was chosen arbitrarily, we have 

E[?7] < -(l + o(l)) as 6 ^ oo 
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This also motivates the need for conditions on finiteness of higher order moments of to obtain 
upper bound on the moments of the detection delay. 



From the above theorem, the following corollary easily follows. 



Corollary 1 ( l\T2y): If the likelihood ratios are such that (26 1, (27 1 and (28 1 are satisfied for some 



r > 1, then for the Shiryaev stopping time Tj defined in ( |2T] ) 

g a\ 



inf E[(r-r) 



n{Ts 



q + d 

A similar result can be concluded for the conditional moments as well. 



as a — )■ 0, for all m < r. 



(34) 



C. Performance analysis for i.i.d. model with geometric prior 
We now present the second order asymptotic analysis of the Shiryaev algorithm for the i.i.d. model, 



provided in |12| using the tools from nonlinear renewal theory introduced in Section II-C 



When the observations {Xn} are i.i.d. conditioned on the change point, condition (26l is satisfied and 

^ k+t 

- ^L>(/i||/o) a.s. P„ whent^oo Vn, 



where Z)(/i||/o) is the K-L divergence between the densities /i and /o (see Definition |4]). From Thereom 3.2 
it follows that for Oq, = 

PFA(rs) < a. 



Also, it is shown in | [T2| that if 

0<Z?(/i||/o) <oo and < Z)(/o||/i) < oo, 
then conditions (27i and (28 1 are satisfied, and hence from Corollary [T] 

inf E[(r - r)+]™ ~ E[(rs - r)+^'" ' 



as a 



0. 



(35) 



,^(/i||/o) + |log(l-p)|. 
Note that the above statement provides the asymptotic delay performance of the Shiryaev algorithm. 
However, the bound for PFA can be quite loose and the first order expression for the ADD ( |35] ) may not 
provide good estimate for ADD if the PFA is not very small. In that case it is useful to have a second 
order estimate based on nonlinear renewal theory, as obtained in |12]. 



First note that ( |25] ) for the i.i.d. model will reduce to 

1 + A 



A^ 



n+l 



l-p 



(36) 
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with Aq = 0. Now, let Z„ = log A„, and recall that Yfc = log fl^xl) ■ Then, it can be shown that 



n-l 



Zn = J][yfc + |log(l-p)|]+log(e^«+p)+ J]log(l + e-^'=p) 



k=l 



k=l 



= ^Yk + n\logil-p)\+7]n. (37) 

k=l 

Therefore the Shiryaev algorithm can be equivalently written as 

Ts = inf{n > I ■ Zn > b} . 

We now show how the asymptotic overshoot distribution plays a key role in second order asymptotic 
analysis of PFA and ADD. Since, p^s > ^ implies that > log j^a ~ logo = b, we have. 



1 



> 



1 + e^^T 1 + a 



PFA(Ts) = E[l - pr,] 



Thus, 



E 
E 



1 + e/ 



< E [: 



e 



1 1 

e^^s 1 + e"^^s 



> E 


1 a 


_e^"s 1 + a_ 





E [e^^^s] (1 + o(l)) as a oo. 



PFA(rs) = E[e~^-s](l + o(l)) = e-*E[e-(^-s-^)](l + o(l)) as 6 ^ oo. 

= e"''E[e"(^-s-f')|r >r](l + o(l)) as 6 ^ oo. 

and we see that PFA is a function of the overshoot when Z„ crosses a from below. 
Similarly, 

Zrs = XI ^fc + "^sl log(l - p)\+Vts- 



k=l 



Following the developments in Section II-C if the sequence y„ satisfies some additional conditions, then 
we can write 

Ei[Z,] -Ei[r/,] 



Ei[r] 



|log(l-/5)|+Z)(/i||/o) 
b + Ei[Zr-b]-Ei[r]r] 
\\og{l-p)\+Dif4fo)- 



It is shown in 1 12 1 that r/„ is a slowly changing sequence, and hence the distribution of Zr^ — 6, as 6 



oo, is equal to the asymptotic distribution of the overshoot when the random walk X]fe=i '^k+^l log(l— p)| 

^As explained in jl2) , this analysis is facilitated if we restrict to the worst-case detection delay, which is obtained by 
conditioning on the event that the change happens at time 1. 
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crosses a large positive boundary. We define the following quantities: the asymptotic overshoot distribution 
of a random walk 

R{x) = lim P ( V Yfc + Tsl log(l - p)\-b<x] (38) 



Its mean 



and its Laplace transform at 1 



\k=l 



oo 

K = I xdR{x) 
'o 



C= / e-^dR{x). 
Jo 



(39) 



(40) 



Also, the sequence {Yn} satisfies some additional conditions and hence the following results are true. 
Theorem 3.3 ( [J2j): If {Yn}s are nonarithmetic then rjn is a slowly changing sequence. Then by 
Theorem 12.71 



This implies 



If in addition E[Y'^] < oo then 



lim ¥{Zr, - b<x) = R{x). 

6— ^-oo 



PFA(rs) ~ C e"'' as 6 ^ oo. 
b + K-Ei[rj] 



+ o(l), as 6 — > oo. 



|log(l-p)|+Z)(/i||/o) 
where rj is the a.s. limit of the sequence {rjn}. 

In Table |ll| we compare the asymptotic expressions for PFA and ADD given in Theorem |3.3| with 
simulations. As can be seen in the table, the asymptotic expressions get more accurate as PFA approaches 
0. 



TABLE II: /o ~ 7^(0, 1), /i ~ 1), p = 0.01 





PFA 


ADD 


b 


Simulations 


Analysis 


Simulations 


Analysis 






Theorem 3.3 




Theorem 3.3 


1.386 


1.22 xlO"^ 


1.39 xlO"^ 


6.93 


10.31 


2.197 


5.85 xlO"^ 


6.19 xlO"^ 


8.87 


11.9 


4.595 


5.61 xlO"^ 


5.63 xlO"^' 


13.9 


16.6 


6.906 


5.59 xlO"* 


5.58 xlO^'' 


18.59 


21.13 


11.512 


5.6 xlO"'' 


5.58 xlO-" 


27.64 


30.16 
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45 




4 5 6 7 8 9 10 11 



I log (PFA) I > 

Fig. 3: ADD - PFA trade-off curve for the Shiryaev algorithm: p = 0.01, /o = 7V(0, 1), /i = 7\A(0.75, 1) 



In Fig. [3] we plot the ADD as a function of log(PFA) for Gaussian observations. For a PFA constraint 
of a that is small, h | log(a)|, and 

Ann ~ w 1 ~ |log(«)l 

giving a slope of | iog(i,p)|+D(/,||/„) to the trade-off curve. 

When I log(l — p)\ L'(/i||/o), the observations contain more information about the change than the 
prior, and the tradeoff slope is roughly TJij^]^- Oi^ the other hand, when | log(l — /o)| S> D(/i||/o), the 
prior contains more information about the change than the observations, and the tradeoff slope is roughly 
I iog(i-p)| • latter asymptotic slope is achieved by the stopping time that is based only on the prior: 

r = inf {n > 1 : P(r > n) < a] . 



This is also easy to see from ([14]). With D{fi\\fo) small, L{X) ^ 1, and the recursion for pk reduces to 

Pn+l =Pn + {l- Pn)p, PO = 0. 

Expanding we get p„ = pY12=o(^ ~ P)^- "^^'^ desired expression for the mean delay is obtained from 
the equation pr = I — a. 

IV. MiNIMAX QUICKEST CHANGE DETECTION 

When the distribution of the change point is not known, we may model the change point as a 
deterministic but unknown positive integer 7. A number of heuristic algorithms have been developed 
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in this setting. The earliest work is due to Shewhart |[T|, Q, in which the log likelihood based on the 
current observation is compared with a threshold to make a decision about the change. The motivation 
for such a technique is based on the following fact: if X represents the generic random variable for the 
i.i.d. model with /q and /i as the pre- and post-change p.d.fs, then 

Eoo[logL(X)] = -Z?(/o||/i)<0 and Ei [logL(X)] = Z?(/i||/o) > 0. (41) 

where as defined earlier L{x) = /i(x)//o(x), and Eqo and Ei correspond to expectations when 7 = oo 
and 7 = 1, respectively. Thus, after 7, the log likelihood of the observation X is more likely to be above 
a given threshold. Shewhart's method is widely employed in practice due to its simplicity; however, 
significant performance gain can be achieved by making use of past observations to make the decision 
about the change. Page [3] proposed such an algorithm that uses past observations, which he called the 
CuSum algorithm. The motivation for the CuSum algorithm is also based on ( |4T] ). By the law of large 
numbers, "^l^^^og L{Xi) grows to 00 as n — )• 00. Thus, if Sn = '^i=i^og L{Xi) is the accumulated 
log likelihood sum, then before 7, Sn has a negative drift and evolves towards —00. After 7, 5„ has a 
positive drift and climbs towards 00. Therefore, intuition suggests the following algorithm that detects 
this change in drift, i.e.. 



where b > 0. Note that 



Tc = inf <; n > 1 : ( 5„ - min Sk] >b}. (42) 

l<fc<n 



Sn — min Sh = max > log,L(Xi) = max > logL(Xj). 

l<k<n 0<k<n ^ l<k<n+l 

i=k+l i=k 

Thus, Tc can be equivalently defined as follows. 
Algorithm 2 (CuSum algorithm): 

Tc = inf{n > 1 : W„ > 6}. (43) 

where 

n 

Wn= max y log L{Xi) (44) 

i=k 

The statistic Wn has the convenient recursion: 

Wn+l = {Wn + logL(X„+i)) + , Wo = 0. (45) 
It is this cumulative sum recursion that led Page to call Tq the CuSum algorithm. 
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The summation on the right hand side (RHS) of ( [44| ) is assumed to take the value when k = n + 1. 
It turns out that one can get an algorithm that is equivalent to the above CuSum algorithm by removing 
the term A; = n + 1 in the maximization on the RHS of (|44l), to get the statistic: 



Cn= max VlogL(X,) (46) 

i=k 

The statistic C„ also has a convenient recursion: 

Cn+l = (Cn)+ + logL(X„+i), Co = 0, (47) 

Note that unlike the statistic Wn, the statistic C„ can be negative. Nevertheless it is easy to see that both 
Wn and Cn cross will cross a positive threshold b at the same time (sample path wise) and hence the 
CuSum algorithm can be equivalently defined in terms of C„ as: 

Tc = inf {n > 1 : Cn > 6} . (48) 

An alternative way to derive the CuSum algorithm is through the maximum likelihood approach i.e., to 
compare the likelihood of {F < n} against {F > n}. Formally, 

- - inf / n M ■ "^^^i^fe^" foW Utk MX,) \ 



Cancelling terms and taking log in ( [49| ) gives us ([48]) with b = log B. 
See Fig. |4] for a typical evolution of the CuSum algorithm. 

Although, the CuSum algorithm was developed heuristically by Page [3], it was later shown in ||6|, |[7|, 
(HJ and pO) , that it has very strong optimality properties. In this section, we will study the CuSum and 
related algorithms from a fundamental viewpoint, and discuss how each of these algorithms is provably 
optimal with respect to a meaningful and useful optimization criterion. 

Without a prior on the change point, a reasonable measure of false alarms is the mean time to false 
alarm, or its reciprocal, which is the false alarm rate (FAR): 

FAR(r) = -\- (50) 

Finding a uniformly powerful test that minimizes the delay over all possible values of 7 subject to a FAR 
constraint is generally not possible. Therefore it is more appropriate to study the quickest change detection 
problem in a minimax setting in this case. There are two important minimax problem formulations, one 
due to Lorden |[6| and the other due to PoUak ||9|. 
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Fig. 4: Typical Evolution of the CuSum algorithm. Threshold 6 = 8 and change point 7 = 60. 



In Lorden's formulation, the objective is to minimize the supremum of the average delay conditioned 
on the worst possible realizations, subject to a constraint on the false alarm rate. In particular, if we 
defin^] 

WADD(r) = sup ess sup E„ [(r - n)+\Xi, . . . . (51) 

n>l 

and denote the set of stopping times that meet a constraint a on the FAR by 

Va = {T: FAR(r) < q} (52) 
we have the following problem formulation due to Lorden. 

Lorden's Problem: For a given a, find a stopping time r G Va to minimize WADD(r). (53) 



For the i.i.d. setting, Lorden showed that the CuSum algorithm (43 1 is asymptotically optimal for 



Lorden's formulation ( [53l ) of as a — )• 0. It was later shown in |[7| and Q that the CuSum algorithm is 



actually exactly optimal for ( |53) . Although the CuSum algorithm enjoys such a strong optimality property 
under Lorden's formulation, it can be argued that WADD is a somewhat pessimistic measure of delay. 
A less pessimistic way to measure the delay was suggested by PoUak ||9j|: 

CADD(r) =sup E„[r-n|r > n]. (54) 

n>l 

^Lorden defined WADD with (r — n + 1)"*" inside the expectation, i.e., he assumed a delay penalty of 1 if the algorithm stops 
at the change point. We drop this additional penalty in our definition in order to be consistent with the other delay definitions 
in this chapter. 
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for all stopping times r for which the expectation is well-defined. 
Lemma 1: 

WADD(r) > CADD(r). 
Proof: Due to the fact that r is a stopping time on 

{T>n] = {t <n-lY e ^2, . . . , Xn-i) 

Therefore, for each n 

ess sup E„ [(r — n)^\Xi, . . . > En,[(r — n)^\T > n] = E„[r — n\T > n] 

and the lemma follows. ■ 
We now state PoUak's formulation of the problem that uses CADD as the measure of delay. 

Pollak's Problem: For a given a, find a stopping time r G to minimize CADD(t). (55) 



PoUak's formulation has been studied in the i.i.d. setting in [9] and |20|, where it is shown that algorithms 
based on the Shiryaev-Roberts statistic (to be defined later) are within a constant of the best possible 
performance over the class Da, as a — )• 0. 

Lai fTO} studied both ( [53] ) and ( [55] ) in a non-i.i.d. setting and developed a general minimax asymptotic 
theory for these problems. In particular, Lai obtained a lower bound on CADD(r), and hence also on 
the WADD(r), for every stopping time in the class D^, and showed that an extension of the CuSum 
algorithm for the non-i.i.d. setting achieves this lower bound asymptotically as a — )• 0. 



In Section IV-A we introduce a number of alternatives to the CuSum algorithm for minimax quickest 



change detection in the i.i.d. setting that are based on the Bayesian Shiryaev algorithm. We then discuss 



the optimality properties of these algorithms in Section [TV-B While we do not discuss the exact optimality 



of the CuSum algorithm from Q or |[8|, we briefly discuss the asymptotic optimality result from |6|. 
We also note that the asymptotic optimality of the CuSum algorithm for both (|53]) and dSSll follows from 



the results in the non-i.i.d. setting of pO|, which are summarized in Section IV-C 



A. Minimax Algorithm's Based on the Shiryaev Algorithm 



Recall that the Shiryaev algorithm is given by (see ( [20| ) and ([23])): 

Ts = inf {n > 1 : Rn^p > a} 

where 

^ k=l i=k 
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Also recall that Rn^p has the recursion: 



Rn+l,p — — ; —L{Xn+l), Ro,p — 0. 

1 — p 

Setting /9 = in the expression for Rn^p we get the Shiryaev-Roberts (SR) statistic pT 



Rn = Y,llL{Xi) (56) 

k=l i=k 

with the recursion: 

Rn+i = (1 + Rn)LiXn+i), Ro = 0. (57) 
Algorithm 3 (Shiryaev-Roberts (SR) Algorithm): 

T,^ = mf{n>l:Rn>B, Rq = 0} . (58) 
It is shown in (^\ that the SR algorithm is the limit of a sequence of Bayes tests, and in that limit 



it is asymptotically Bayes efficient. Also, it is shown in |20| that the SR algorithm is second order 
asymptotically optimal for ( |55| ), as a — 0, i.e., its delay is within a constant of the best possible delay 
over the class Va- Further, in |[22j, the SR algorithm is shown to be exactly optimal with respect to a 
number of other interesting criteria. 

It is also shown in |[9| that a modified version of the SR algorithm, called the Shiryaev-Roberts-PoUak 



(SRP) algorithm, is third order asymptotically optimal for ( [55| ), i.e., its delay is within a constant of the 
best possible delay over the class Va, and the constant goes to zero as a — )• 0. To introduce the SRP 
algorithm, let be the quasi-stationary distribution of the SR statistic i?„ above: 

F{Ro <x) = Q^{x) = lim Po(^n. < x\t,^ > n). 

n— >oo 

The new recursion, called the Shiryaev-Roberts-PoUak (SRP) recursion, is given by, 

Rn+l = (1 + Rn)L{Xn+l) 

with Rq is distributed according to . 
Algorithm 4 (Shiryaev-Roberts-Pollak (SRP) Algorithm): 

T,^, = m{{n>l:R^>B }. (59) 



Although the SRP algorithm is strongly asymptotically optimal for PoUak's formulation of ( pS] ), in 
practice, it is difficult to compute the quasi-stationary distribution Q^. A numerical framework for 
computing efficiently is provided in [23 1. Interestingly, the following modification of the SR algorithm 
with a specifically designed starting point Rq = r > Ois found to outperform the SRP procedure uniformly 
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over all possible values of the change point |20|. This modification, referred to as the SR-r procedure, 
has the recursion: 



■n+l 



[1 + Rn)L{Xn+i), Ro 



Algorithm 5 (Shiryaev-Roberts-r (SR-r) Algorithm): 

T,,_^ = mi{n>l:Rl>B}. (60) 



It is shown in |20| that the SR-r algorithm is also third order asymptotically optimal for ( |55] ), i.e., its 
delay is within a constant of the best possible delay over the class Va, and the constant goes to zero as 
a — )■ 0. 



Note that for an arbitrary stopping time, computing the CADD metric ([54]) involves taking supremum 
over all possible change times, and computing the WADD metric ( [ST] ) involves another supremum over 
all possible past realizations of observations. While we can analyze the performance of the proposed 
algorithms through bounds and asymptotic approximations (as we will see in Section IV-B[ it is not 
obvious how one might evaluate CADD and WADD for a given algorithm in computer simulations. This 
is in contrast with the Bayesian setting, where ADD (see (|7])) can easily be evaluated in simulations by 
averaging over realizations of change point random variable T. 



Fortunately, for the CuSum algorithm ( [43] ) and for the Shiryaev-Roberts algorithm ( 158] ), both CADD 
and WADD are easy to evaluate in simulations due to the following lemma. 
Lemma 2: 

CADD(tc) = WADD(rc) = Ei \(t^ - 1)] . (61) 

CADD(ts,) = WADD(ts,) = El [(ts, - 1)] . (62) 

Proof: The CuSum statistic Wn (see ( [44| )) has initial value and remains non-negative for all n. 
If the change were to happen at some time n > 1, then the pre-change statistic Wn-\ is greater than 
or equal 0, which equals the pre-change statistic if the change happens at n = 1. Therefore, the delay 
for the CuSum statistic to cross a positive threshold h is largest when the change happens at n = 1, 
irrespective of the realizations of the observations, Xi,X2, . . . Therefore 

WADD(tc) = sup ess sup E„ [(tc - n)+|Xi, . . . ,X„_i] = Ei [(r^ - 1)+] = Ei \(t^ - 1)] . 

n>l 

and 

CADD(rc) = sup E„[t - n|rc > n] = Ei [(tc - l)|rc > 1] = Ei \{t^ - 1)] . 

n>l 



October 23, 2012 



DRAFT 



28 



This proves ( [6T] ). A similar argument can be used to establish ([62]). ■ 
Note that the above proof crucially depended on the fact that both the CuSum algorithm and the Shiryaev- 
Roberts algorithm start with the initial value of 0. Thus it is not difficult to see that Lemma |2] does not 
hold for the SR-r algorithm, unless of course r = 0. Lemma [2] holds partially for the SRP algorithm 
since the initial distribution makes the statistic stationary in n. As a result E„[tsrp — n|rsRp > n\ 
is the same for every n. However, as mentioned previously, is difficult to compute in practice, and 
this makes the evaluation of CADD and WADD in simulations somewhat challenging. 

B. Optimality Properties of the Minimax Algorithms 

In this section we first show that the algorithms based on the Shiryaev-Roberts statistics, SR, SRP, and 
SR-r are asymptotically optimal for PoUak's formulation of ( [55] ). We need an important theorem that is 
proved in [[22|. 



Theorem 4.1 ( /|22j/j.- If the threshold in the SR algorithm ( [58] ) can be selected to meet the constraint 
a on FAR with equality, then 

• E^=l^n[(T-n)+] 

= argmm — — 

Proof: We give a sketch of the proof here. Note that 

oo oo 

J]E„[(r-n)+] = J^^'^t^-^l^^^lIP'ooCr >n) 
71=1 n=l 

oo 

< CADD(r)^Poo(r > n) 

n=l 

= CADD(r)Eoo[r], 

and hence is well defined for all stopping times for which CADD and FAR exist. The first part of the 
proof is to show that 

oo oo 
J];E„[(rsR - n)+] = min J]E„[(r - n)+]. 

71=1 71=1 

This follows from the following result of ||9|. If 

Jx,p{r) = min (e [(r - r)+] + A P(t < T] 



with r having the geometric distribution of (111 with parameter p, and 



T"A,p = arg min JA,p(r). (63) 
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Then for a given Tsr (with a given threshold), there exists a sequence {{Xi, pi)} and with Aj — )• A* and 
Pi ^ such that r^^^p^ converge to Tsr, as i — )• oo. Thus, the SR algorithm is the limit of a sequence of 
Bayes tests. Moreover, 

limsup -^'P^ _ I 

p^O,A-s>A* 1 — J\,p{Tsr) 

By ([63|), for any stopping time r, it holds that 



Now by taking the limit p — 0, A — A* on both sides, using the fact that for any stopping time r ||22 

1 - ^A,p(r) 



A*Eoo[r] - ^E[(r-n)+] as p ^ 



n=l 

and using the hypothesis of the theorem that the FAR constraint can be met with equality by using Tsr, 
we have the desired result. 

The next step in the proof is to show that it is enough to consider stopping times in the class Va that 
meet the constraint of a with quality. The result then follows easily from the fact that is optimal with 
respect to the numerator m ^^-^ — ■ ■ 

We now state the optimality proof for the procedures SR, SR-r and SRP. We only provide an outline 
of the proof to illustrate the fundamental ideas behind the result. 

Theorem 4.2 ( [20]): If IEi[log < oo, and log j^^l^"! is nonarithmetic then 

1) 

inf CADD(r) > , +k + o(l) as q ^ (64) 



where k is a constant that can be characterized using renewal theory |[T8|. 

2) For the choice of threshold B = ^, FAR(rsR) < a, and 

I lo£ Oil 

CADD(rsR) = '„,,' + C + 0(1) as a ^ 
= inf CADD(r) + 0(1) as a ^ 

where ( is again a constant that can be characterized using renewal theory fisl. 

3) There exists a choice for the threshold B = B^ such that FAR(tsrp) < a(l + o(l)) and 

I loE C(\ 

CADD(rsRp) = '/,,„' , +k + o(l) as a ^ 
= inf CADD(r) + o(l) as a ^ 0. 



(65) 



(66) 
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4) There exists a choice for the threshold B = such that FAR(tsrp) < q(1 + o(l)) and a choice 
for the initial point r such that 

logal „ , , 

— + K + o(l) as a — 

(67) 



CADDfr., 



DihWh) 



= inf CADD(t) +o(l) as a ^ 0. 

Proof: To prove that the above mentioned choice of thresholds meets the FAR constraint, we note 
that {EI^ — n — r} is a martingale. Specifically, {Rn — n} is a martingale and the conditions of theorem 
2.3| are satisfied ||24|. Hence, 

r,J = 0. 



Since, EooI-Rtsr] > B, for B 



we have 
FAR(rs,) 



1 



1 

< 

- B 



a. 



For a description of how to set the thresholds for Tsr_^ and Tsrp, we refer the reader to |20|. 

The proofs of the delay expressions for all the algorithms have a common theme. The first part of 

^,E„[(7--n)+] ^ ^^^^^ ^^^^^ CADD(r). 

E?1iSUP„IE„[t - n|r > n]W^{T > j) 



these proofs is based on Theorem 



4.1 



We first show that 



CADD(r) = supE„[r - n|r > n] 



> 



From Theorem 



4.1 



inf CADD(r) > 



ET=i^j[r-j\r>j]F^{T>j) 

Er=i%[(r-i)+] 

inimizii 

E^=ilEn[(TsR -n)- 



since Tsr is optimal with respect to minimizing ^"'^^"^[^J] ^ we have 



The next step is to use nonlinear renewal theory (see Section pI-C[ ) to obtain a second order approximation 

ve did 

log a I 



for the right hand side of the above equation, as we did for the Shiryaev procedure in Section III-C 



n) 



+ k + o(l) as a 



0. 



KooM i^(/i||/o) 
The final step is to show that the CADD for the SR-r and SRP procedures are within o(l), and the 
CADD for SR procedure is within 0(1) of this lower bound ( [64] ). This is done by obtaining second 
order approximations using nonlinear renewal theory for the CADD of SR, SRP, SR-r procedures, and 



get (65 1, (66 1 and ([67|), respectively. 
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J L V°° E [(t— 

It is shown in |22| that jg" r i — — is also equivalent to the average delay when the change happens 



at a "far horizon": 7 — )■ oo. Thus, the SR algorithm is also optimal with respect to that criterion. 

The following corollary follows easily from the above two theorems. Recall that in the minimax setting, 
an algorithm is called third order asymptotically optimum if its delay is within an o(l) term of the best 
possible, as the FAR goes to zero. An algorithm is called second order asymptotically optimum if its 
delay is within an 0(1) term of the best possible, as the FAR goes to zero. And an algorithm is called 
first order asymptotically optimal if the ratio of its delay with the best possible goes to 1, as the FAR 
goes to zero. 



Corollary 2: Under the conditions of Theorem 4.1 the SR-r (60 1 and the SRP (p9b algorithms are 



third order asymptotically optimum, and the SR algorithm is second order asymptotically optimum for 

to 



the PoUak's formulation ( pS] ). Furthermore, using the choice of thresholds specified in Theorem 4.1 
meet the FAR constraint of a, the asymptotic performance of all three algorithms are equal up to first 
order: 

CADD(ts,) ~ CADD(ts,p) ~ CADD(rs,_,0 



log a I 



DihWfo) 



Furthermore, by Lemma |2] we also have: 

I log a\ 



WADDfr., 



DihWfoY 

In Q the asymptotic optimality of the CuSum algorithm ([43]) as a — )■ is established for Lorden's 



formulation of (53 1. First, ergodic theory is used to show that choosing the threshold b = \ loga| ensures 



FAR(tc) < a. For the above choice of threshold B = \ loga|, it is shown that 

WADD(rc) < ' (1 + 0(1)) as a ^ 0. 

Then the exact optimality of the SPRT |[25| is used to find a lower bound on the WADD of the class Pq,: 



inf WADD(t) > ' (1 + 0(1)) as a ^ 0. 

These arguments establish the first order asymptotic optimality of the CuSum algorithm for Lorden's 



formulation. Futhermore, as we will see later in Theorem |4.3| Lai 1 10 1 showed that: 



inf WADD(r) > inf CADD(r) > f 1 + o(r 

Thus by Lemma [2j the first order asymptotic optimality of the CuSum algorithm extends to Pollak's 
formulation ( [55] ), and we have the following result. 
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f^=N(0,l) f^=N(0.75,l) 




7 7.5 8 8.5 9 9.5 10 

-log (FAR) > 



Fig. 5: ADD - PFA trade-off curve for the CuSum algorithm: /o = M{0, 1), /i = AA(0.75, 1) 



Corollary 3: The CuSum algorithm ( |43] ) with threshold 6 = |loga| is first order asymptotically 
optimum for both Lorden's and PoUak's formulations. Furthermore, 

CADD(Te) = WADD(re) ~ 

In Fig. |5} we plot the trade-off curve for the CuSum algorithm, i.e., plot CADD as a function of 
-log FAR. Note that the curve has a slope of |/o). 

C. General Asymptotic Minimax Theory 



In pO], the non-i.i.d. setting earlier discussed in Section III-B is considered, and asymptotic lower 



bounds on the WADD and CADD for stopping times in Va are obtained under quite general conditions. 



It is then shown that the an extension of the CuSum algorithm ( [43| ) to this non-i.i.d. setting achieves this 
lower bound asymptotically as a — )• 0. 

Recall the non-i.i.d. model from Section |III-B| When the process evolves in the pre-change regime, 
the conditional density of Xn given is fQ^n{Xn\^i^'^)- After the change happens, the conditional 

density of X„ given X^~^ is given by Also 

h,{x,\x{-^) 

The CuSum algorithm can be generalized to the non-i.i.d. setting as follows: 
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Algorithm 6 (Generalized CuSum algorithm): Let 

Cn = max y^li- 

l<k<n ^ 
i=k 

Then the stopping time for the generaUzed CuSum is 



Tg = inf {n > 1 : C„ > b} (68) 



Then the following result is proved in 1 10|. 

Theorem 4.3: If there exists a positive constant I such that the {li}s satisfy the following condition 



n+rn 



lim sup ess sup P„ max N Yi > 7(1 + (5)n 



then, as a — 



=0 V5>0 (69) 



inf WADD(r) > inf CADD(r) > ( 1 + o(l) ) . (70) 



Further 

Koo[ro] > e^ 

and under the additional condition of 

^lim sup ess sup ( t^^ Yi > I — 6 



m+t 

1 =0 V(5>0, (71) 



m>n 



we have 



CADD(rG) < WADD(to) < y ( 1 + o(l) ) as 6 ^ oo. 



Thus, if we set 6 = | loga|, then 



and 



WADD(r,) < + 



which is asymptotically equal to the lower bound in ( [70| ) up to first order. Thus is first-order asymp- 
totically optimum within the class Va of tests that meet the FAR constraint of a. 

Proof: We only provide a sketch of the proof for the lower bound since it also based on the idea of 
using Chebyshev's inequality. The fundamental idea of the proof is to use Chebyshev's inequality to get 
a lower bound on any arbitrary stopping time r from Va, such that the lower bound is not a function of 
r. The lower bound obtained is then a lower bound on the CADD for the entire family D^. 
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Let La and Va be positive constants. To show that 

sup E„[t — n|r > n] > La{l + o{l)) as a — )• 

n>l 

it is enough to show that there exists n such that 

E,„[r - n|r > n] > Lq,(1 + o(l)) as a 0. 

This n is obtained from the following condition. Let m be any positive integer. Then if r G Va, there 
exists n such that 

Foo{t > n) > and Poo(r < n + m|r > fc) < ma. (72) 



We use the n that satisfies the condition ^T2\ . Then, by Chebyshev's inequality 

^n{T -n> La\T > n) < (L„)"^E„[r - n|r > n]. 

We can then write 

E„[r — n|r > n] > Lq, Pn(''" — n > La\T > n). 

Now it has to be shown that Fn{T — n > La\T > n) ^ I uniformly over the family Va- To show this, 
we condition on Va- 



(T-n<La\T>n) = Pn \ T - n < La] y^^Yj < Va\T > 



n 



+Fnir-n<La;Y,Y,> Va\T > 



n 



The trick now is to use the hypothesis of the theorem and choose proper values of Va and La such that 
the two terms on the right hand side of the above equations are bounded by terms that go to zero and 
are not a function of the stopping time r. We can write 



i T ~ n < La ; > I T > n j < ess sup P„ I max 



In flOl , it is shown that if we choose 

log a I 



La = {l- 6) 

and 



Va = il-6^)\\oga\ 
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then the condition (69 1 ensures that the above probabiUty goes to zero. The other term goes to zero by 



using a change of measure argument and using condition ( |72| ): 

P„ I T - n< La;^Yi < V^jr > n | < e^°Poo(r <n + La\T > n). 



V. Relationship Between The Models 



We have discussed the Bayesian formulation of the quickest change detection problem in Section III 



and the minimax formulations of the problem in Section IV The choice between the Bayesian and the 
minimax formulations is obvious, and is governed by the knowledge of the distribution of T. However, it 
is not obvious which of the two minimax formulations should be chosen for a given application. As noted 
earlier, the formulations of Lorden and PoUak are equivalent for low FAR constraints, but differences 



arise for moderate values of FAR constraints. Recent work by Moustakides |26| has connected these 
three formulations and found possible interpretations for each of them. We summarize the result below. 

Consider a model in which the change point is dependent on the stochastic process. That is, the 
probability that change happens at time n depends on {Xi, . . . , X„}. Let 

TTn = n^ = n\Xi,...,Xn). 

The distribution of F thus belongs to a family of distributions. Assume that while we are trying to find 
a suitable stopping time r to minimize delay, an adversary is searching for a process {vr„} such that the 
delay for any stopping time is maximized. That is the adversary is trying to solve 

J(r) = sup E[t-T\t > T]. 
M 

It is shown in p6| that if the adversary has access to the observation sequence, and uses this information 



to choose TTn, then J(r) becomes the delay expression in Lorden's formulation ( [53| ) for a given r, i.e.. 



J(t) = sup ess sup E„,[(t — n)^\Xi, . . . , Xn^i]. 

n>l 

However, if we assume that the adversary does not have access to the observations, but only has access 



to the test performance, then J(r) is equal to the delay in PoUak's formulation ( [55] ), i.e.. 



J(r) = sup E„[t — n|r > n]. 

n>l 



Finally, the delay for the Shiryaev formulation ( [TO] ) corresponds to the case when 7r„ is restricted to only 
one possibility, the distribution of F. 
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VI. Variants and Generalizations of the Quickest Change Detection Problem 

In the previous sections we reviewed the state-of-the-art for quickest change detection in a single 
sequence of random variables, under the assumption of complete knowledge of the pre- and post-change 
distributions. In this section we review three important variants and extensions of the classical quickest 
change detection problem, where significant progress has been made. We discuss other variants of the 
change detection problem as part of our future research section. 

A. Quickest Change Detection with unknown pre- or post-change distributions 

In the previous sections we discussed quickest change detection problem when both the pre- and 
post-change distributions are completely specified. Although this assumption is a bit restrictive, it helped 
in obtaining recursive algorithms with strong optimality properties in the i.i.d. setting, and allowed the 
development of asymptotic optimality theory in a very general non-i.i.d. setting. In this section, we 
provide a review of the existing results for the case where this assumption on the knowledge of the pre- 
and post-change distributions is relaxed. 

Often it is reasonable to assume that the pre-change distribution is known. This is the case when 
changes occur rarely and/or the decision maker has opportunities to estimate the pre-change distribution 
parameters before the change occurs. When the post-change distribution is unknown but the pre-change 
distribution is completely specified, the post-change uncertainty may be modeled by assuming that the 
post-change distribution belongs to a parametric family {Pe}- Approaches for designing algorithms in 
this setting include the generalized likelihood ratio (GLR) based approach or the mixture based approach 
In the GLR based approach, at any time step, all the past observations are used to first obtain a maximum 
likelihood estimate of the post-change parameter 9, and then the post-change distribution corresponding 
to this estimate of 6 is used to compute the CuSum, Shiryaev-Roberts or related statistics. In the mixture 
based approach, a prior is assumed on the space of parameters, and the statistics (e.g., CuSum or Shiryaev- 
Roberts), computed as a function of the post change parameter, are then integrated over the assumed 
prior. 

In the i.i.d setting this problem is studied in \^ for the case when the post-change p.d.f. belongs to 
a single parameter exponential family {fe}, and the following generalized likelihood ratio (GLR) based 
algorithm is proposed: 




(73) 
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If ea = I jpg^i , and Cq, is of the order of | log a\, then it is shown in |6| that as the FAR constraint a — )• 0, 



FAR(ro) <a(l + o(l)), 
and for each post-change parameter 6* / 0, 

WADD^(r,) 



DifeJoY 

Here, the notation WADD^[-] is used to denote the WADD when the post-change observations have p.d.f 
fg. Note that £^j^'f^-^ is the best one can do asymptotically for a given FAR constraint of a, even when 
the exact post-change distribution is known. Thus, this GLR scheme is uniformly asymptotically optimal 



with respect to the Lorden's criterion ( [53) . 

In |j27J, the post-change distribution is assumed to belong to an exponential family of p.d.f. s as above, 
but the GLR based test is replaced by a mixture based test. Specifically, a prior G(-) is assumed on the 
range of 6 and the following test is proposed: 



inf < n > 1 : max 

l<fc<n 



>-}■ (74) 



For the above choice of threshold, it is shown that 

FAR(rM) < a, 

and for each post-change parameter ^ / 0, if is in a set over which G{-) has a positive density, as the 
FAR constraint a — ;> 0, 

WADOnrJ ~ ^ 

^[Je,Jo) 

Thus, even this mixture based test is uniformly asymptotically optimal with respect to the Lorden's 
criterion. 

Although the GLR and mixture based tests discussed above are efficient, they do not have an equivalent 
recursive implementation, as the CuSum test or the Shiryaev-Roberts tests have when the post-change 
distribution is known. As a result, to implement the GLR or mixture based tests, we need to use the 
entire past information {Xi,--- ,Xn), which grows with n. In ||Tol, asymptotically optimal sliding- 
window based GLR and mixtures tests are proposed, that only used a finite number of past observations. 
The window size has to be chosen carefully and is a function of the FAR. Moreover, the results here 
are valid even for the non-i.i.d. setting discussed earlier in Section |IV-C[ Recall the non-i.i.d. model 
from Section III-B with the prechange conditional p.d.f. given by fo^n{Xn\X]^~^), and the post-change 
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conditional p.d.f. given by /i ^] 
algorithms are given, respectively, by: 



To = inf < n > 1 : max sup 

n—mo,<k<n—m' 



Then the generalized CuSum ( [68] ) GLR and mixture based 

fB,i{Xi\X[ ^) 



fo,i{Xi\Xl 



and 



Tm = inf < n > 1 : max 

I n—m„<k<n 



(75) 



(76) 



It is shown that for a suitable choice of Ca, rua and m^, under some conditions on the likelihood ratios 
and on the distribution G, both of these tests satisfy the FAR constraint of a. In particular, 

FAR(fM) < a, 

and as a — )• 0, 

FAR(fo) <a(l + o(l)). 
Moreover, under some conditions on the post-change parameter 6, 

WADD^[fJ ~ WADD^[f, 



log a I 



DifeJo)' 

Thus, under the conditions stated, the above window-limited tests are also asymptotically optimal. We 
remark that, although finite window based tests are useful for the i.i.d. setting here, we still need to store 
the entire history of observations in the non-i.i.d. setting to compute the conditional densities, unless the 



dependence is finite order Markov. See [28] and |10| for details. 

To detect a change in mean of a sequence of Gaussian observations in an i.i.d. setting, i.e., when 



/o ^ AA(0, 1) and /i ~ M{fi, 1), ^ / 0, the GLR rule discussed above (75 1 (with nia = n and = 1) 
reduces to 



inf < n > 1 : max 

0<fc<n 



t y/{n-k) 



> b 



(77) 



This test is studied in |29| and performance of the test, i.e., expressions for FAR(rt,) and Ei[t„] are 
obtained. In 1 28 1, it is shown that a window-limited modification of the above scheme ( [77] ) can also be 
designed. 

When both pre- and post-change distributions are not known, again GLR based or mixture based tests 



have been studied and asymptotic properties characterized. In |30|, this problem has been studied in the 
i.i.d setting (Bayesian and non-Bayesian), when both pre- and post-change distributions belong to an 
exponential family. For the Bayesian setting, the change point is assumed to be geometric and there are 
priors (again from an exponential family) on both the pre- and post-change parameters. GLR and mixture 
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based tests are proposed that have asymptotic optimaUty properties. For a survey of other algorithms when 
both the pre- and post-change distributions are not known, see | [28| . 

In |3T|, rather than developing procedures that are uniformly asymptotically optimal for each possible 
post-change distribution, robust tests are characterized when the pre- and post-change distributions belong 
to some uncertainty classes. The objective is to find a stopping time that satisfies a given false alarm 
constraint (probabihty of false alarm or the FAR depending on the formulation) for each possible 
pre-change distribution, and minimizes the detection delay (Shiryaev, Lorden or the PoUak version) 
supremized over each possible pair of pre- and post-change distributions. It is shown that under some 
conditions on the uncertainty classes, the least favorable distribution from each class can be obtained, 
and the robust test is the classical test designed according to the lease favorable distribution pair. It is 



shown in |31| that although robust test is not uniformly asymptotically optimal, it can be significantly 
better than the GLR based test of |6] for some parameter ranges and for moderate values of FAR. The 
robust solution also has the added advantage that it can be implemented in a simple recursive manner, 
while the GLR test does not admit a recursive solution in general, and may require the solution to a 
complex nonconvex optimization problem at every time instant. 

B. Data-efficient quickest change detection 

In the classical problem of quickest change detection that we discussed in Section [lll| a change in 
the distribution of a sequence of random variables has to be detected as soon as possible, subject to a 
constraint on the probability of false alarm. Based on the past information, the decision maker has to 
decide whether to stop and declare change or to continue acquiring more information. In many engineering 
applications of quickest change detection there is a cost associated with acquiring information or taking 
observations, e.g., cost associated with taking measurements, or cost of batteries in sensor networks, etc. 



(see |32| for a detailed motivation and related references). In |32|, Shiryaev's formulation (Section [III]) is 
extended by including an additional constraint on the cost of observations used in the detection process. 
The observation cost is captured through the average number of observations used before the change 
point, with the understanding that the cost of observations after the change point is included in the delay 
metric. 

In order to minimize the average number of observations used before F, at each time instant, a decision 
is made on whether to use the observation in the next time step, based on all the available information. 
Let Sn G {0, 1}, with S'„ = 1 if it is been decided to take the observation at time n, i.e. Xn is available 
for decision making, and 5'„ = otherwise. Thus, Sn is an on-off (binary) control input based on the 



October 23, 2012 



DRAFT 



40 



information available up to time k — 1, i.e., 

Sn = ^J'k-l{Ik-l), A: = 1,2,... 
with fi denoting the control law and / defined as: 



^1; • • • ; ^1 ; • • • ; 



(S) 

Here, represents Xi if Si = 1, otherwise Xi is absent from the information vector In- 

Let ip = {r, /uo, • • . , /Ur-i} represent a policy for data-efficient quickest change detection, where r 



is a stopping time on the information sequence {In}- The objective in |32| is to solve the following 
optimization problem: 



minimize ADD(V') = E [(r - T) 



■4, 

subject to PFA(V') = P(r < T) < a, 

"mm(T,r-l) 

ANO(^) = E 



(78) 



and 



k=l 



</3. 



Here, ADD, PFA and ANO stand for average detection delay, probability of false alarm and average 
number of observations used, respectively, and a and /3 are given constraints. 
Define, 

p„ = P(r<A:|/„). 



Then, the two-threshold algorithm from p2[ is: 

Algorithm 7 (DE-Shiryaev: 'ip{A,B)): Start with pQ = and use the following control, with B < A, 
for k> 0: 



S, 



k+l 



ifpn<B 

1 ifpn>B 

T = inf {k > 1 : pn > A} . 
The probability pn is updated using the following recursions: 

Pn = + (1 - Pn)p if Sk+1 = 



(79) 



Pk+1 



p"„L(Xfc+i) 



if Sk+i = 1 



With L(X,,+i) = /i(X,,+i)//o(Xfe+i). 

With B = the DE-Shiryaev algorithm reduces to the Shiryaev algorithm. When the DE-Shiryaev 



algorithm is employed, the a posteriori probability p„ typically evolves as depicted in Fig. 6a It is shown 
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11 



-log{PFA) »■ 

(a) Evolution of p„ for /o ~ A/'(0, 1), /i ~ A/'(0.5, 1), and (b) Trade-off curves comparing performance of tiie DE- 
p — 0.01, with thresholds A = 0.65 and B — 0.2. Shiryaev algorithm with the Fractional Sampling scheme 

when B is chosen to achieve ANO equal to 50 % of mean 
time to change. Also /o ~ A/'(0, 1), /i ~ A/'(0.75, 1), and 
p = 0.01. 

Fig. 6: Evolution and performance of the DE-Shiryaev algorithm. 



in p2l that for a fixed /3, as a — ;> 0: 



ADDWAi^))~ „(^^.j;7;°'l(,_^)| as.^0 (80, 



and 



PFA('0(^, B)) ~ « (^^ e-'^dRix)^ as a ^ 0. (81) 

Here, R{x) is the asymptotic overshoot distribution of the random walk X]fc=i [1°S ^(^fe) + I log(l — 
p)\], when it crosses a large positive boundary. It is shown in |12| that these are also the performance 
expressions for the Shiryaev algorithm. Thus, the PFA and ADD of the DE-Shiryaev algorithm approach 
that of the Shiryaev algorithm as a — 0, i.e., the DE-Shiryaev algorithm is asymptotically optimal. 

The DE-Shiryaev algorithm is also shown to have good delay-observation cost trade-off curves: for 
moderate values of probability of false alarm, for Gaussian observations, the delay of the algorithm 
is within 10% of the Shiryaev delay even when the observation cost is reduced by more than 50%. 
Furthermore, the DE-Shiryaev algorithm is substantially better than the standard approach of fractional 
sampling scheme, where the Shiryaev algorithm is used and where the observations to be skipped are 



determined a priori in order to meet the observation constraint. (See Fig. 6b ) 

In most practical applications, prior information about the distribution of the change point is not 
available. As a result, the Bayesian solution is not directly applicable. For the classical quickest change 
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detection problem, an algorithm for the non-Bayesian setting was obtained by taking the geometric 



parameter of the prior on the change point to zero (see Section IV-Ai. Such a technique cannot be 
used in the data-efficient setting. This is because when an observation is skipped in the DE-Shiryaev 
algorithm in |32|, the a posteriori probability is updated using the geometric prior. In the absence of 
prior information about the distribution of the change point, it is by no means obvious what the right 
substitute for the prior is. A useful way to capture the cost of observations in a minimax setting is also 
needed. 

In 1 33 1, the minimax formulation of fO"] is used to propose a minimax formulation for data-efficient 
quickest change detection. We observe that in the two-threshold algorithm ^(A, i?), when the change 
occurs at a far horizon, it is the fraction of observations taken before change that is controlled. This 
insight is used to formulate a duty cycle based metric to capture the cost of taking observations before 
the change point. Also, we note that the duration for which observations are not taken in the algorithm 



in 1 32 1, is also a function of the undershoot of the a posteriori probability when it goes below the 
threshold B. Given the fact that for the DE-Shiryaev algorithm, has the interpretation of the the 
likelihood ratio of the hypotheses "ifi : F < n" and "i^o : T > n", the undershoots essentially carry the 



information on the likelihood ratio. It is shown in |33 1 that this insight can be used to design a good test 
in the non-Bayesian setting. This algorithm is called the DE-CuSum algorithm and it is shown that it 
inherits good properties of the DE-Shiryaev algorithm. The DE-CuSum algorithm is also asymptotically 
optimal in a sense similar to ( fSO] ) and ( |8T] ), has good trade-off curves, and performs significantly better 
than the standard approach of fractional sampling. 



C. Distributed sensor systems 

In the previous sections, we provided a summary of existing work on quickest change detection 
and classification in the single sensor (equivalently, centralized multisensor) setting. For the problem of 
detecting biological and chemical agents, the setting that is more relevant is one where there is a set of 
distributed sensors collecting the data relevant for detection, as shown in Fig. [7] Based on the observations 
that the sensors receive, they send messages (which could be local decisions, but not necessarily) to a 
fusion center where a final decision about the hypothesis or change is made. 

Since the information available for detection is distributed across the sensors in the network, these 
detection problems fall under the umbrella of distributed (or decentralized) detection, except in the 
impractical setting where all the information available at the sensors is immediately available without 
any errors at the fusion center. Such decentralized decision making problems are extremely difficult. 
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Channel 



Fusion Center 



Change 
Decision 



Fig. 7: Change detection using distributed sensors 



Without certain conditional independence assumptions across sensors, the problem of finding the optimal 
solutions, even in the simplest case of static binary hypothesis testing, is computationally intractable 134]- 



|39|. Decentralized dynamic decision making problems, of which the quickest change detection problem 
is a special case, are even more challenging due to the fact that information pattern in these problems is 
non-classical, i.e., the different decision makers have access to different pieces of information pO]. 
The problem of decentralized quickest change detection in distributed sensor systems was introduced 



in 1 41 1, and is described as follows. Consider the distributed multisensor system with L sensors, com- 
municating with a fusion center shown in Fig. [v] At time n, an observation is made at sensor 5^. 
The changes in the statistics of the sequences {X^} are governed by the event. Based on the information 
available at time n, a message is sent from sensor Se to the fusion center. The fusion center may 
possibly feedback some control signals to the sensors based on all the messages it has received so far. 
For example, the fusion center might inform the sensors how to update their local decision thresholds. 
The final decision about the change is made at the fusion center. 

There are two main approaches to generating the messages at the sensors. In the first approach, the 
sensors can be thought of simply quantizing their observations, i.e., U!^ is simply a quantized version of 
X^. The goal then is to choose these quantizers over time and across the sensors, along with a decision 
rule at the fusion center, to provide the best tradeoff between detection delay and false alarms. In the 
second approach, the sensors perform local change detection, using all of their observations, and the 
fusion center combines these decisions to make the final decision about the change. 
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The simplest observation model for the decentralized setting is one where the sensors are affected by 
the change at the same time, and the observations are i.i.d. in time at each sensor and independent across 



sensors in both the pre-change and post-change modes. This model was introduced in |41 1, and studied in 



a Bayesian setting with a geometric prior on the change point for the case of quantization at the sensors. 



It was shown that, unlike the centralized problem for which the Shiryaev test is optimal (see Section III i. 



in the decentralized setting the optimization problem is intractable in even for this simple observation 



model. Some progress can be made if we allow for feedback from the fusion center |41 1. Useful results 
can be obtained in the asymptotic setting where the probability of false (Bayesian formulation) or false 
alarm rate (minimax formulation) go to zero. These results can be summarized as follows (see, e.g., [12| , 



|42|, gSI for more details): 

• It is asymptotically optimum for the sensors to use stationary monotone likelihood ratio quantizers, 
i.e., the sensors use the same quantization functions at all times, and the quantization regions 
are obtained by dividing the likelihood ratio of the observations into intervals and assigning the 
quantization levels to them in increasing order. 

• The optimum quantization thresholds at the sensors are chosen to maximize the K-L divergence 
between the post-change and pre-change distributions at the sensors. 

• For fixed stationary quantizers at the sensors, the fusion center is faced with a centralized quickest 
change detection problem. Therefore, depending on the optimization criterion (Bayes or minimax), 
asymptotically optimum change detection procedures can be designed using the techniques described 
in Sections [nil and |TV] 

• The tradeoff between delay and false alarms is governed by the K-L divergence of the quantized 
observations at the output of the sensors, and hence the first order asymptotic performance with 
quantization is at best equal to that without quantization. 

For the case where the sensors make local decisions about the change, it is reasonable to assume that 
the local detection procedures use (asymptotically) optimum centralized (single sensor) statistics. For 
example, in the Bayesian setting, the Shiryaev statistic described in Algorithm [T] can be used, and in the 



minimax setting one of the statistics described in Section IV-A can be used depending on the specific 
minimax criterion used. The more interesting aspect of the decision-making here is the fusion rule used 
to combine the individual sensor decisions. There are three main basic options that can be considered 



T^m'- the fusion center stops and declares the change as soon as one of the sensors' statistics crosses 



October 23, 2012 



DRAFT 



45 



its decision threshold. 

• Tman- the sensors stop taking observations once their local statistics cross their thresholds, and the 
fusion center stops and declares the change after all sensors have stopped. 

• Tali: the sensors continue taking observations even if their local statistics cross their thresholds, and 
the fusion center stops and declares the change after all the sensor statistics are above their local 
thresholds simultaneously. 

It was first shown by |[44| that the t^u procedure using CuSum statistics at the sensors is globally first 



order asymptotically optimal under Lorden's criterion (53 1 if the sensor thresholds are chose appropriately. 



That is, the first order performance is the same as that of the centralized procedure that has access to all 



the sensor observations. A more detailed analysis of minimax setting was carried out in |45|, in which 



procedures based on using CuSum and Shiryaev-Roberts statistics at the sensors were studied under 



PoUak's criterion ( [55] ). It was again shown that Tan is globally first order asymptotically optimal, and that 
Tmax and Tmin are not. 

For the Bayesian case, if the sensors use Shiryaev statistics, then both r^ax and Tan can be shown to 



be globally first order asymptotically optimal, with an appropriate choice of sensor thresholds |43|, |46|. 
The procedure Tmin does not share this asymptotic optimality property. 

Interestingly, while tests based on local decision making at the sensors can be shown to have the same 
first order performance as that of the centralized test, simulations reveal that these asymptotics "kick in" 
at unreasonably low values of false alarm probability (rate). In particular, even schemes based on binary 
quantization at the sensors can perform better than the globally asymptotically optimum local decision 
based tests at moderate values of false alarm probability (rate) |]43|. These results point to the need for 
further research on designing procedures that perform local detection at the sensors that provide good 
performance at moderate levels of false alarms. 

D. Variants of quickest change detection problem for distributed sensor systems 

In Section VI-C[ it is assumed for the decentralized quickest change detection problem, that the change 



affects all the sensors in the system simultaneously. In many practical systems it is reasonable to assume 
that the change will be seen by only a subset of the sensors. This problem can be modeled as quickest 
change detection with unknown post-change distribution, with a finite number of possibilities. A GLRT 
based approached can of course be used, in which multiple CuSum tests are run in parallel, corresponding 
to each possible post-change hypotheses. But this can be quite expensive from an implementation view 
point. In [47], a CuSum based algorithm is proposed in which, at each sensor a CuSum test is employed. 
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the CuSum statistic is transmitted from the sensors to the fusion center, and at the fusion center, the 
CuSum statistics from all the sensors are added and compared with a threshold. This test has much lower 
computational complexity as compared to the GLR based test, and is shown to be asymptotically as good 
as the centralized CuSum test, as the FAR goes to zero. Although this test is asymptotically optimal, the 
noise from the sensors not affected by change can degrade the performance for moderate values of false 
alarm rate. In |48|, this work of pTl , is extended to the case where information is transmitted from the 
sensors only when the CuSum statistic at each sensor is above a certain threshold. It is shown that this 



has the surprising effect of suppressing the unwanted noise and improving the performance. In |49|, it 
is proposed that a soft-thresholding function should be used to suppress these noise terms, and a GLRT 
based algorithm is proposed to detect presence of a stationary intruder (with unknown position) in a 
sensor network with Gaussian observations. A similar formulation is described in ||50l. 

The Bayesian decentralized quickest change detection problem under an additional constraint on the 
cost of observations used is studied in [51]. The cost of observations is captured through the average 
number of observations used until the stopping time and it is shown that a threshold test similar to the 
Shiryaev algorithm is optimal. Recently, this problem has been studied in a minimax setting in ||52| and 



asymptotically minimax algorithms have been proposed. Also, see |53| and |54| for other interesting 



energy-efficient algorithms for quickest change detection in sensor networks. 

VII. Applications of quickest change detection 

As mentioned in the introduction, the problem of quickest change detection has a variety of applications. 
A complete list of references to applications of quickest change detection can be quite overwhelming 
and therefore we only provide representative references from some of the areas. For a detailed list of 
references to application in areas such as climate modeling, econometrics, environment and public health. 



finance, image analysis, navigation, remote sensing, etc., see [13] and [55 1. 

1) Statistical Process Control (SPC): As discussed in the introduction, algorithms are required that 
can detect a sudden fault arising in an industrial process or a production process. In recent years 
algorithms for SPC with sampling rate and sampling size control have also been developed to 
minimize the cost associated with sampling [ [56} , [ |57| . See [ [58} , and [ [59} and the references therein 
for some recent contributions. 

2) Sensor networks: As discussed in [|32j, quickest change detection algorithms can be employed in 



sensor networks for infrastructure monitoring [60|, or for habitat monitoring |61 1. Note that in these 
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applications, the sensors are deployed for a long time, and the change may occur rarely. Therefore, 
data-efficient quickest change detection algorithms are needed (see Section VI-B| ). 



3) Computer network security: Algorithms for quickest change detection have been applied in the 



detection of abnormal behavior in computer networks due to security breaches |62|, |63|, |64|. 
4) Cognitive radio: Algorithms based on the CuSum algorithm or other quickest change detection 
algorithms can be developed for cooperative spectrum sensing in cognitive radio networks to detect 
activity of a primary user. See |53[ , | |65| , and ||67|. 



5) Neuroscience: The evolution of the Shiryaev algorithm is found to be similar to the dynamics of 



the Leaky Integrate-and-Fire model for neurons 1 68 1 . 



6) Social Networks: It is suggested in |69| and |70| that the algorithms from the change detection 



literature can be employed to detect the onset of the outbreak of a disease, or the effect of a 
bioterroist attack, by monitoring drug sales at retail stores. 

VIII. Conclusions and future directions 

In this article we reviewed the state-of-the-art in the theory of quickest change detection. We saw that 
while exactly or nearly optimal algorithms are available only for the i.i.d. model and for the detection of a 
change in a single sequence of observations, asymptotically optimal algorithms can be obtained in a much 
broader setting. We discussed the uniform asymptotic optimality of GLR and mixture based tests, when 
the post-change distribution is not known. We discussed algorithms for data-efficient quickest change 
detection, and showed that they are also asymptotically equivalent to their classical counterparts. For the 
decentralized quickest change detection model, we discussed various algorithms that are asymptotically 
optimal. We also reviewed the asymptotic optimality theory in the Bayesian as well as in the minimax 
setting for a general non-i.i.d. model, and showed that extensions of the Shiryaev algorithm and the 
CuSum algorithm to the non-i.i.d. setting are asymptotically optimal. Nevertheless, the list of topics 
discussed in this article is far from exhaustive. 

Below we enumerate possible future directions in which the quickest change detection problem can 
be explored. We also provide references to some recent articles in which some research on these topics 
has been initiated. 

1) Transient change detection: It is assumed throughput this article that the change is persistent, i.e., 
once the change occurs, the system stays in the post-change state forever. In many applications it 
might be more appropriate to model change as transient, i.e., the system only stays in the post- 



change state for a finite duration and then returns to the pre-change state; see e.g., |71 1, |72| and 
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|73|. In this setting, in addition to false alarm and delay metrics, it may be of interest to consider 
metrics related to the detection of the change while the system is still in the change state. 
2) Change propagation: In applications with multiple sensors, unlike the model assumed in Sec- 



tion 



VI-C it may happen that the change does not affect all the sensors simultaneously |74|. The 
change may propagate from one sensor to the next, with the statistics of the propagation process 
being known before hand | [75| . 
3) Multiple change points in networks: In some monitoring application there may be multiple change 
points that affect different sensors in a network, and the goal is to exploit the relationship between 



the change points and the sensors affected to detect the changes |76|. 
4) Quickest change detection with social learning: In the classical quickest change detection problem, 
to compute the statistic at each time step, the decision maker has access to the entire past observa- 
tions. An interesting variant of the problem is quickest change detection with social learning, when 
the time index is replaced by an agent index, i.e., when the statistic is updated over agents and not 
over time, and the agents do not have access to the entire past history of observations but only to 



some processed version (e.g., binary decisions) from the previous agent; see | [77| , | |78| and | |79[ . 
5) Change Detection with Simultaneous Classification: In many applications, the post-change distri- 
bution is not uniquely specified and may come from one of multiple hypotheses Hi, . . . ,Hm, in 
which case along with detecting the change it of interest to identify which hypothesis is true. See, 



e.g., 1 80 1 and |81|. 



6) Synchronization issues: If a quickest change detection algorithm is implemented in a sensor network 
where sensors communicate with the fusion center using a MAC protocol, the fusion center might 
receive asynchronous information from the sensors due to networking delays. It is of interest to 
develop algorithms that can detect changes while handling MAC layer issues; see, e.g., (50}. 
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